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Abstract 

This paper is related to our previous works [1] [2] on the error estimate 
of the averaging technique, for systems with one fast angular variable. In the 
cited references, a general method (of mixed analytical and numerical type) 
has been introduced to obtain precise, fully quantitative estimates on the 
averaging error. Here, this procedure is applied to the motion of a satellite 
in a polar orbit around an oblate planet, retaining only the J2 term in the 
multipole expansion of the gravitational potential. To exemplify the method, 
the averaging errors are estimated for the data corresponding to two Earth 
satellites; for a very large number of orbits, computation of our estimators 
is much less expensive than the direct numerical solution of the equations of 
motion. 
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1 Introduction. 



Celestial mechanics often requires approximation techniques to compute the motion 
of bodies over very long times. Giving reliable error estimates on these techniques 
is a very practical problem, of not simple solution. 

In this paper we apply to a typical astronomical problem the general scheme of 
[1] [2] to estimate the error of the averaging principle, in the case of one fast angular 
variable d (in the one-dimensional torus) and many slow variables / = (P) (the 
actions). The application we consider refers to the motion of a satellite around an 
oblate planet: of course, oblateness gives a perturbation of the Keplerian, purely 
radial gravitational potential. There is a classical way to express the perturbed 
potential as a series of spherical harmonics, where the multipole moments of the 
planetary mass distribution appear as coefficients. In the case of axial symmetry, 
this series only involves zonal harmonics and the £-th term contains a multipole co- 
efficient Je (for a reminder, see subsection 3E and references therein). For a slightly 
oblate planet, all these coefficients are small if they are defined in a convenient, 
dimensionless fashion; often, J\ ~ and one retains only the 1 = 1 term of the 
previous expansion; the corresponding equations of motion for the satellite are the 
classical " J 2 problem". 

In this paper the J 2 problem is considered in the restricted case of polar orbits; 
this problem involves a fast angular variable 6 (the angle between the polar axis and 
the radius vector of the satellite) and three slow variables I = (P, E, Y) representing 
the parameter, the eccentricity and the argument of the pericenter for the osculating 
Keplerian ellipse. In principle, all these variables are unknown functions of the 
"physical" time t; however, it is customary to consider as a "time" variable the 
angle itself or, rather, the ratio t of the angle to 2n. With this choice of the time 
variable, by definition the angle evolves with the law 0(t) = [2nt} (where [ ] indicates 
equivalence mod.27r); for obvious reasons, t will be called the orbit counter. The 
evolution of I = (P, E, Y) is described by a set of equations of the form 



where e is a small parameter proportional to J 2 (see again subsection 3E). Indepen- 
dently of this specific application, a system of evolution equations like (1.1), with 
any number of unknown functions I = (l l )i=i,...,d and a small parameter e > 0, will 
be called in the sequel a perturbed periodic system. To any such system one can 
associate an averaged system, replacing the functions /' in (1.1) with their averages 
over the angle 6 = [27rt]; the solution J = (J 1 ) of the latter is in a function of r := et. 

The difference between the solutions I of (1.1) and J of the averaged system, 
on a time scale of order 1/e, can be estimated in a fully quantitative way with the 
general method of [1] [2]; this is what we are going to do in the present work, for 
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the periodic system describing the polar J 2 problem. As in the cited papers, the 
expression "fully quantitative" means that our final estimate will have the form 

|r(t) - J*(et)| < en l (ei) for t e [0, U/e) , (1.2) 

where n l : [0, U) — > [0, +oo) are computable functions, and U is a constant, whose 
choice determines quantitatively the interval of observation; we note that U/e is the 
total number of orbits. 

The availability of an algorithm to construct the estimators n* is the main differ- 
ence between the approach of [2] and the classical, qualitative result I*(t) — J l (et) = 
0(e), see, e.g., [3] [4] [5]. We are aware of the existence of higher order versions 
of the averaging method, in which the error behaves like 0(e p ) for t G [0, 0(l/e)) 
(p = 2, 3, see again [4] [5]). However, the classical theory of higher order averag- 
ing gives little more than the qualitative error estimate 0(e p ); perhaps an extension 
of our methods could give quantitative estimators of the form e p (et) , but this 
problem is left to future work. 

Let us return to the functions n* in Eq. (1.2). Our method to compute them for 
the J2 problem is mainly analytical, but, in its final steps, it requires the numerical 
solution of an ODE on the interval [0, U); however, this numerical computation is 
much faster than the direct numerical solution of (1.1) (a fact already appearing in 
other applications, different from the J2 problem, considered in our previous papers). 

The general treatment proposed in this paper for the polar J2 problem is sub- 
sequently exemplified, choosing for e and for the initial conditions Jo = (Pq, Eq, Yq) 
the values corresponding to the Earth and to the Polar and Cos-B satellites. In 
this case, the estimators n* have been computed up to 60000 orbits (i.e., one or two 
centuries), spending few seconds of CPU time on a PC. To test their reliability, we 
have computed by direct numerical integration the differences I l (t) — J l (et); this is 
a more expensive operation, that we have been able to perform up to 3000 orbits 
with the same PC; our estimators are quite satisfactory on this interval, since the 
functions 1 1— > en l (et) are close to the envelopes of the rapidly oscillating functions 

Of course, the treatment of real satellites of the Earth should include, besides 
the J2 gravitational term, other minor perturbations such as: gravitational forces 
corresponding to higher order moments of the Earth, atmospheric dragging, solar 
wind, tidal effects due to the Moon gravity. All these effects could be treated, 
giving rise again to a system of the form (1.1) with more complicated perturbation 
components ef l (even for non polar orbits: in general, if the analysis is not restricted 
to the orbits in a fixed plane, the actions are not three but five). Presumably, the 
corresponding averaged system and our estimators n* for its error could be computed 
as well; this is left to future work. In spite of the limitation to the polar J 2 problem, 
we think that the results of this paper have some interest, because they show that 
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a general method for error estimates works over very long times on a non-trivial 
problem. 

To conclude, we describe in few words the organization of the paper. In Section 
2, we show how to compute the error estimators for the averaging of any periodic 
system (1.1); this illustration specializes to the periodic case the slightly more general 
setting of [1] [2] for one-frequency systems. In Section 3, the basic facts on the 
motion of a satellite in a fixed plane II are reviewed; the outcomes are the equations 
of motion for (I*) = (P, E, Y), for any perturbation of the Keplerian potential keeping 
the satellite on II. In Section 4, we specialize the previous equations to the polar J 2 
problem; the averaged system is solved, and we construct the algorithm to compute 
the error estimators n* (i = p, e, y). In Section 5, we perform the computations for 
the Polar and Cos-B satellites of the Earth. Some details on the constructions of 
Sections 3-4 are presented in Appendices A-D. 



2 Averaging of periodic systems. 

2A. Introducing the problem. We consider an open set A of (the space of the 
actions) H d and the one-dimensional torus T: 

A = {/ = (/'), ; 4 C R d , T := R/2ttZ = {0} ; (2.1) 

the natural projection of the real axis on the torus will be written [ ] : R — > T, 
x i — ^ [x]. Now, let 

/ = !./•' ), , ,/ G C m (A x T, R d ), (J, 0) » /(/, 0) (2.2) 

(with m ^ 2); we fix 

J GA, £G(0,+oo), (2.3) 
and consider the Cauchy problem 

f { = s/(l, [27rt]) , 1(0) = J : (2.4) 

the maximal solution (in the future) is a C m+1 function I : [0, t max ) C R — > A, 
t I— > l(t) (with t max G (0, +oo]). Here and in the sequel, typeface symbols are 
employed for functions of t (or r = et)); this allows, for example, to distinguish the 
above mentioned function I from a point I of A. 
The averaged system associated to (2.4) is 

% = 7( J) , J(0) = Jo , (2.5) 
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f = (n =1 ,..., d eC m (A,R d ) , 



I ^ 



f(T):=^J dfif(I,-ff) ; (2.6) 



the maximal solution (in the future) is a C m+1 function J : [0, r max ) cRhA,th 
J(t). Throughout the paper, we are interested in binding the difference I(t) — J(et). 

2B. Connections with the framework of [1] [2]. In these papers, we have 
discussed the Cauchy problem 



with / as in (2.2) and g G C m (A x T,R), to G C m (A,R), u(I) ^ for all / G A, 
I G A, -&Q G T, e G (0, +oo), the maximal solution being (I, 0) : [0, t max ) — > A x T. 
Precise quantitative estimates have been derived on the difference I(t) — J(et), where 
J is the solution of (2.5). 

Clearly, the perturbed periodic system (2.4) is equivalent to a particular case of 
(2.7). In fact, let us choose 

:= 2vr , g(I, #) := for all / G A, # G T ; tf := ; (2.8) 

then, the equation for 6 in (2.7) is fulfilled with 



and the equation for I in (2.7) becomes, with this position, the Cauchy problem 
(2.4). This remark allows to apply the general results of [1] [2] to the problem (2.4); 
in the sequel, we report from these papers the minimal elements enabling one to read 
independently the present paper, i.e.: (i) some basic assumptions required by our 
approach, and the definitions of certain auxiliary functions; (ii) the final proposition 
from [2], estimating l(t) — J(et) in terms of these auxiliary functions. 

2C. A first set of assumptions and auxiliary functions. We use the tensor 
spaces and the operations on tensors already employed in our previous works; in 
particular, the spaces T} (R d ), T 2 (R d ) and T^(R d ) are made, respectively, by the 
the (1, 1), (2, 0) and (1, 2) tensors over R d ; in the three cases, these are represented 
as families of real numbers srf = (^), = (e^ 1 - 7 ), ^ = (^ fc ) k — 1, d). 

From here to the end of the section the function /, the initial datum I and the 
parameter e of Eqs. (2.2) (2.3) are fixed. We stipulate the following: 

(1) s,v G C m (A x T,R d ), p,q,w G C m -\A x T,R d ), u G C m - 2 (A x T, ~R d ) and 
^# G C m ~ 2 (A, T}(R d )) are the auxiliary functions uniquely defined by the equations 




o 



(2.7) 



9(1) := M 



(2.9) 




(2.10) 



s = 




v(I, 0) = for all I G A; 



(2.11) 
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»>:=£/; ("2) 

<9u> 

P = P + 2tt— , w(J,0) = for all Je A; (2.13) 
aw 

<9w <9 2 7- /<97\ 2 

" : =w /; ■"■■=w'-{<n) ■ (214) 

The barred symbols /, s, p denote the averages of /, s,p over the angle, in the sense 
of (2.6); d/dl and d 2 /dl 2 indicate the Jacobian and the Hessian with respect to the 
variables I = (P). Explicit formulas for s,v and w are 

s = z-z, *(/,*):= ^ j^dff (/(/,<?')-/(/)); (2-15) 



16) 



«(/, i?) := ^ s(J, tf') ; «;(/, 0) := ^ jf (W (p(J, tf') - p(/)) . (2. 

The above definitions of s,p, ...,*/# specialize the general prescriptions of [1] [2] to 
the case (2.8), of interest in this paper. Similar comments could be added in the 
sequel, but they will not be repeated. 

(2) We introduce the open set 

A t := {(/, 51) G A x R d | [/, / + 51} C A} , (2.17) 

where [/, I + 51} is the closed segment in R d with the indicated extremes. From 
now on, <S G C m ~ 2 (Af, T}(R d )) and G C m - 2 (A t , T*(R d )) are two functions such 
that, for all {1,51) G A t , 

p(I + 5I) = p(I)+&(I,6I)6I , (2.18) 
— — (97 1 

f(I + 5I) = f(I) + ^ I (I)5I+-^(I,5I)5I 2 , je* k (I,6I) = Jt£ j (I,6I). (2.19) 

and J4? are not uniquely determined, if d > 1: possible choices, given by Taylor's 
formula, are 

&(I,5I):= [ dx^(I+x5I), Jt?(I,6I) :=2 / . (2.20) 

(3) From now on, [0, [/) (£/ G (0, +oo]) is a fixed interval where the solution J of 
the averaged system (2.5) is assumed to exist. We denote with R : [0, U) — > T}(R d ), 
t i— > R(t) and K : [0, U) — > R d , r i— > K(t) the solutions of the Cauchy problems 

f = f(J)R, R(0) = 1„; (2.21) 



dK d f , , 

^ = ^y(J)K + P(J) , K(0) = 0; 



(2.22) 



these are C m , and exist on the whole interval [0, U) due to the linearity of the 
above differential equations. By standard arguments, one proves that R(r) is an 
invertible matrix for all r G [0, U), and derives for K the explicit formula K(r) = 

2D. Some more auxiliary functions. We maintain the assumptions and nota- 
tions of the previous items (1) (2) (3). In [2] we have introduced a second set of 
auxiliary functions, related to s,p, ....R, K and to any separating system of seminorms 
on H d . For simplicity, here we consider the seminorms giving the absolute values of 
each component of vectors and tensors in R d ; what follows is an adaptation of [2] to 
this choice. From now on, k range in {1, d} and we use Einstein's summation 
convention on repeated indices. 

(4) For J G R d and g = (g*) G [0, +oo} d , we put 



B(J, g) := {I G R d | \ f - f \ < ^ V i} . 
(5) p = (p l ) e C([0, U), [0, +oo} d ) is a function such that 

B(J(t),p(t)) C A for tG [0,17). 

We put 

T p := {(r, r) G [0, U) x [0, +oo) d | r* < p l (r) V % } 



^2.23) 



(2.24) 



(2.25) 



(6) a*, 6* in C 2 (r p , [0, +oo)) and c*,d},e} fe = e 1 ^ G C 1 ^, [0, +oo)) are functions 
such that for any r G [0, [/), 57 G 77(0, p(r)), # G T, 



s( J(r) + 5 J, t?) - R(r)s(/ , #„) " K(r) 



(w(J(t) + 5 J, 0) - ||-(J(t)) „( J(r) + *J, 0))' 



(«(J(r) + <5J, 0) - ^(J(r))(w + g)(j(r) + 5 J, 0) 

-^(J(t))u(J(t)+ 5.7,0))] ^c\t,\5J\) , 

|^(j(r),5J)K4(r,|5J|) , 

|^l(j(r),5J)K4 fe (r,|5J|)- 
In the above, one always intends 

M := (\SJ\%=i,..., d • 



^a*(r,M) , 



<b\r,\5J\) , 



(2.26) 

(2.27) 
(2.28) 

(2.29) 
(2.30) 

(2.31) 
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The functions c\ d^, e* fc are assumed to be non-decreasing with respect to the variable 
r, i.e., 

(T,r),(T,r')eT p , r^r' J Vj c l (r,r) ^ c l (r,r') (2.32) 

and similarly for dp and e*- fe . 

(7) Given a*, — , e* fc , we define the functions 



a* G C 2 (r p , [0, +oo)), a l (r, r) := a*(r, r) + eft^r, r) 
fe^^xlO.+oof, [0,+oo)), 



(2.33) 
(2.34) 



1 . 



7 (J, r, 



\r,r) + d)(r,r)f + -e) k (r,r)f£ k . 



(8) i2j G C\[0,U), [0,+oo)) and P] G C([0, £/), [0, +oo)) are functions such that, 
forrG[0,C/), 

|R}(r)Kit!}(r) , KOJ^K P](r) . (2.35) 

2E. The main result on general periodic systems. In [2], we have shown 
how to obtain estimators for I(t) — J(^t) solving a system of integral equations, 
or of equivalent differential equations. Here we only report the final differential 
reformulation, that will be used in the present paper. Keeping in mind the previous 
items (l)-(8), we have the following statement. 

2.1 Proposition, (i) Assume there are = G [0 + oo) d ; (A l j) G [0, +oo) d2 
and a = (a 1 ) G (0 + oo) d , such that 

£ := lb : ,/>: - + a 1 ] C n i =i,..., d (0,p i (0)/e) 



and 



A 



: = max A) < 1/e , 



da 1 



(0,e£) 



^A) fori.,, \....d. feE, 



a*(0,e^) -f*| +eA)a % < a 1 for i = l,...,d. 



Then, there is a unique £ 



such that 



£ GE, a (o,ee )=e . 



(2.36) 

(2.37) 

(2.38) 
(2.39) 

(2.40) 



(ii) With £ as above, let m = (m l ),n = (n l ) G C 1 ([0, U), R d ) solve the Cauchy 
problem 

^ = p;y(., £ „,n) , m*(0) = 0, (2.41) 



n l (0)=£* 

/or all i, with the domain conditions 

< en* < p l , det (l - e-£ (•, en)) > (2.43) 

(m £/ie above, 1 — eda/dr stands for the matrix (5* — eda l /dr^) (i,j = l,...,d), and 
Eq. (2.42) contains the matrix elements of its inverse. We note that (2.4-1) implies 
m % ^ 0). Then, the solution I of the periodic system (2.4) exists on [0, U/e) and 

\l\i)- f{ei)\^e\C{et) for i = 1, d, t e [0, U/e). (2.44) 

2.2 Remarks. For future use, it is worthy to repeat some considerations of [2]. 

(a) The previous Proposition mentions £ , the unique fixed point of the map a(0, e-) 
in the set £ of Eq. (2.36). The proof given in [2] indicates that a(0, e-) : S — > S has 
Lipschitz constant £\A < 1 in the maximum component norm \\z\\ := maxj \z l \, with 
j4. as in (2.37). So, from the standard theory of contractions, we have the iterative 
construction 

£ = lim l n , li any point of E , l n :— a(0, el n -\) for n — 2, 3, ... . (2.45) 

n— »+oo 

For each n ^ 2, 

IKo-ZnlKMr 1 ^^. (2.46) 

One can use the above characterization of £ to compute it numerically; in this case 
one finds the iterates up to a large order n, and then approximates £ with l n . 

(b) In typical cases, and also in the forthcoming application to satellites, the Cauchy 
problem (2.41) (2.42) for the unknown functions n\ m* is solved numerically (by some 
standard package for ODEs). 

2F. A final comment on the functions a 1 , ...,e*. fc . Let us start with a remark 
on b\ c\ dp e l j k . These functions are always multiplied by a small factor e in the main 
statements on \¥{t) — J l (ei)|, i.e., Proposition 2.1; for this reason, in applications it 
is generally sufficient to determine b\ ...,e* fc majorizing roughly the left-hand sides 
ofEqs. (2.27)-(2.30). 

The situation is different for the functions a 1 , which are not multiplied by e; in 
this case, it is important to determine them estimating accurately the left-hand side 
of Eq. (2.26) or, at least, the part of order zero in 8 J. To this purpose, we note the 
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existence of a function G C™" 1 ^ x T, T}(R )) such that, for all (J, 51) G A t 
and G T, 

s(J + 5/, ■&) = s(J, t?) + J^(J, (5/, ; (2.47) 
a solution of this equation is given by Taylor's formula, i.e., 

&(I,5I,-&):= J dx—(I + x5I,$). (2.48) 

To go on, we write 

s( J(r) + <fj, t?) - R(r)s(/ , 7? ) - K(r) (2.49) 

= (s( J(r), #) - R(r)s(/ , # ) - K(r)) + J(r), 5 J, $)5J , 

which decomposes the function in the left-hand side of (2.26) into a zero order 
part in 5J plus a reminder controlled by . Now, assume there are functions 
G C 2 ([0,£7), [0,+oo)) and a) G C 2 (r p , [0, +oo)) such that, for all r G [0,17), 



57 G 5(0, p(r)) and G T, 



s(J(r),7?)-R(r)s(/ ,^o)-K(r) 



^ «(0)( r ) , 



|^(j(r),5J,^)Ka}(r,|5J|) . 
Then, Eq. (2.26) is fulfilled by the function 

a\r,r) := aj 0) (r) + a)(r,r)r j ; 



(2.50) 



(2.51) 



(2.52) 



of course, this definition must be inserted in Eq. (2.33) for 7. In applications one 
will determine al -, binding very accurately the left-hand side of Eq. (2.50), and a*- 
binding more roughly the left-hand side of Eq. (2.51); this is convenient, since in 
Eq. (2.52) for a 1 the terms a* are multiplied by r j and this will finally take the small 
values r J = enP(r). 

The above strategy to determine a\ e* fe will be employed in the application of 
the next sections, concerning satellite motions. 



3 Basic facts on satellite motions. 

3 A. A general setting. A satellite of mass m, position P and velocity P is 
assumed to move around a planet of mass M. As a first approximation, we assume 
the mass of the planet to be uniformly distributed within a sphere of center O; so, 
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the gravitational potential energy (per unit mass) of the planet has the Keplerian 
form 

W = (3-D 

where G ~ 6.674- 10~ n m 3 Kg~ 1 sec -2 is the gravitational constant. The correspond- 
ing gravitational force — mgradV/<-(P) will be referred to as the Keplerian force, and 
all the other forces are assumed to be small perturbations of it. 

The non-Keplerian force acting on the satellite (of gravitational, or any other 
nature) will be written as emf(P, P), where e > is a small dimensionless param- 
eter, arising naturally from the analysis of this perturbation; the parameter and the 
satellite mass are factored out from the perturbation, for convenience. So, the total 
force acting on the satellite is 

F(P, P)=m (-gradV^(P) + ef(P, P)) . (3.2) 

We assume the perturbation to keep the satellite on a fixed plane II passing through 
O, where we introduce the polar coordinates p(P) := \P — 0\, d(P) := angle (in T) 
between a fixed unit vector k and P — O. The equations for the satellite motions in 
this plane can be written as a first-order system for p, $ and their derivatives p, d, 
in the following way: 

%-pP + ™=eQ p (p,#,p,#), | = p, (3.3) 

p 2 — + 2pp& = eQ*(p,0,p,0) , — = 0, (3.4) 

where Q p := f-(dP/dp) and Q$ = f-{dP/dd) are the Lagrangian components of 
the perturbation /. In particular, in the conservative case f(P) = — gradW(P) we 
have 

dW dW 
Q p {p,$) = —fy(p,#) , Q*(P,#) = —Q#(P^) ■ ( 3 - 5 ) 

3B. Kepler elements. Our aim is to re-express the system (3.3) (3.4) after a 
change of coordinates (p, p, i?) i— > (P,E,Y,"&) where P,E,Y are related to the 
unperturbed Kepler problem (as in Eqs. (3.3) (3.4), with Q p = and Q# — 0) . For 
the sake of brevity, let us introduce the abbreviations 

£ = (p,<?,p), = (p,&,p,&) . (3.6) 

For given (£,#), we define 

0£# := unperturbed Kepler orbit in II issuing from the initial datum (£,$) ; (3.7) 
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as usually, we refer to this as the osculating Kepler orbit. Let 

^ '■= I $ > ) is a (nondegenerate) ellipse } 

(circular or rectilinear trajectories being excluded), and define maps 

e ^ ^ P ^ e (o, +oo), E(£, 0) e (0, l), y (£, 0) e T (3.8) 

where P, 77, y are, respectively, the parameter/P, the eccentricity and the argument 
of the pericenter for the ellipse o^; here, 

R := a characteristic length of the problem (3-9) 

(e.g., the planetary radius) . In the above, the standard parameter is divided by R to 
obtain a dimensionless quantity; from now on, P itself will be called the parameter 
and we will refer to (P, E, Y, •&) as the Kepler elements of o^. As well known, S> is 
the set of states with d > and negative unperturbed energy: 

={(£,<?) |0>O, Up 2 + P 2 $ 2 ) - — < 0} ; (3.10) 

2 p 

furthermore, 



y(f , t?) := the solution <¥ E T of the equation (3.13) 
p = l + m ^^_^ such that signsin(tf - <*) = signp . 
To continue, we introduce the notation 

/ = (P, E, Y) = (7 P , I E , I Y ) (3.14) 
and write (£, i?) i— > 7(£, i?) for the correspondence (3.11)-(3.13). The mapping 

^ - (0, +oo) x (0, 1) x T 2 , (£, t?) h- (7(£, i?), t?) (3.15) 
is one-to-one, with inverse 

(0,+oo) x (0,1) x T 2 -> , (3.16) 
(7,0)- , 
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P{I^)--=\{^E^-Y) , (3.17) 



^ (J ' V := V (i + ^ cos (^ - ^)) 2 > ( 3 - 18 ) 

(see Appendix A for more details). For future use, hereafter we give the apocenter 
p+, the pericenter p_ and the orbital period for the unperturbed Kepler motion on 
an ellipse of parameter P and eccentricity E; these are 



RP / R 3 P 3 

Of course, the first two equations can be inverted giving 

P( P+ , P _ ) = ^±£— E(p + ,p.) = ^£. (3.21) 

R{p+ + p-) P+ + P- 

3C. The equations of motions in terms of the Kepler elements. Using the 
transformations (3.11)— (3.19) and Eqs. (3.3) (3.4), one easily obtains a system of 
differential equations for (I, ■&) , namely 



d~F / RP 

3£ + 4 cos(tf - y) + E cos(2tf - 2F) . , 

+£ 27^1F Q » M ' 



+£ E\/GMRP ' 



ill) I GM ,, _, 2 



Q„(/,tf) := (?„«(/,#),#) , (a € {,,, /)}) . (3.26) 

As expected, for Q a = the above equations reproduce the fact that P, E, Y are 
constants of the motion for the unperturbed Kepler problem. 
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From now on, for technical reasons related to the application of the general 
averaging method we will think the "slow" angle Y as an element of R rather than 
T; so, the space of Kepler elements becomes 

A := {/ = (P, E,Y)\Pe (0, +00), E G (0, 1), Y G R} . (3.27) 

Taking this viewpoint amounts to identify the vector field defined by the above 
equations with its lift with respect to the projection (P, E, Y, $) G A x T 1— > 
(P, E, [Y], 1?) G (0, 1) x (0, +00) x T x T (recall that [ ] : R -> T is the quotient 
map). 

3D. From the physical time t to the "orbit counter" t. Let us consider the 
(maximal) solution t G [0, i max ) 1— > (I(t),$(t)) of Eqs. (3.22)-(3.25), with initial data 

7(0) = J = (P , E , Y ) , tf(0) = , (3.28) 

and denote by t( ) : [0, t max ) 1— > R, £ 1— > t(t) the unique C 1 function such that 

t(0) = , #(t) = [27rt(f)] . (3.29) 

Then 

d~& , s ^ dt. . ,„ 
dt W = dt W ' ( } 

and the positivity of dfi/dt ensures the function 1 1— > t(t) to be increasing. The image 
of this function is an interval [0, t max ), and we denote the inverse function by 

t( ) : [0, t max ) -> [0, t max ) , t ^ t(t) . (3.31) 

To go on, let us define 

J(t) := [J(f)]«(t) (3-32) 

for t G [0,t max ). Then di/dt = (dl / dt) / (dt/ dt) = 2ir(dl/dt)/(dti/dt) and Eqs. 
(3.22)-(3.25) with the data (3.28) yield 



dP AneP 2 RQ#(I,ti) 



+7VE P 



dt (1 + Pcos(^-F)) 2 GM 

dE 2neP 2 sm($-Y) R 2 Q P (I,$) 
~dt ~ (1 + Pcos(^-F)) 2 GM 

3P + 4 cos(tf -Y) + E cos(2tf - 2Y) P Q#(1, #) 



P(0) = P , (3.33) 

1?=[27Tt] 



;i + Pcos(^-F)) 2 GM 



(3.34) 
P(0) = P , 



1?=[27Tt] 



dY_ P 2 cos(^-r) p 2 q p (i,i) 

dt ~ P(l + Pcos(^-F)) 2 GM [6 ' 6b) 
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n Psm($-Y)(2 + Ecos($-Y) RQM,$) 
E(l + E cos($ - Y)) 2 GM 



Y(0) = Y ; 

l9=[27Tt] 



we note that t, R 2 Q P (I, d)/{GM) and RQ#(I ',■&) / (GM) are dimensionless quan- 
tities. Eqs. (3.33)-(3.35) are a system of the general form (2.4) for the function 
t i— > 7(t). Eq. (3.29) indicates that t(t) is the number of orbits performed by the 
satellite from time to time t; this explains the name of orbit counter employed for 
the t variable. 

3E. Polar motions of a satellite around an oblate planet. From now on, 
the perturbation is due to a deviation of the planetary mass from the uniform 
distribution inside a sphere. The mass distribution is assumed to be symmetric with 
respect to the polar axis; O is the midpoint between the poles, and we introduce 
a system of Cartesian coordinates x,y,z with origin O, with the z axis (indicated 
by a unit vector k) as polar axis The total gravitational potential P h- > V(P) 
produced by the planet admits a well known expansion in zonal harmonics, namely, 



v(P)- 0X1 



\P-0 

where Pg are the Legendre polynomials, and 



R e ( z(P) \ 

l -Y,\iTTw JtPt \\P=W 

e=i 



(3.36) 



In the above: M is total mass of the planet; R is a typical length expressing the 
size of the planet, usually chosen as the equatorial radius; Pg are the Legendre 
polynomials; dM is the measure describing the mass distribution of the planet. 
We note that Ji, J 2, J3, ■■■ are dimensionless due to the presence of the parameters 
M and R in their definitions; we will refer to them as the (dimensionless) dipole, 
quadrupole, octupole,... coefficients. For a derivation of (3.36) one can refer, e.g., 
to [6]. 

As for the applications of this expansion to (polar or generic) satellite motions, 
there is an enormous literature; we only cite the classical references [5] and [7]. To 
continue, let us recall that 

A(C):=C, ^(C):=^(3C 2 -1); (3.38) 

in general Pe(— £) = (—l) Pg(£), which implies Jg = if I is odd and the mass 
distribution is reflection-invariant with respect to the equatorial plane. Often, J\ ~ 
even without an exact equatorial symmetry; in this case, the first non-negligible 



1 Of course, we are assuming m/M to be so small that O can be regarded at rest. 
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contribution to the expansion (3.36) is the t = 2 term. From now on, we assume 
Ji = and neglect the terms of (3.36) with £ ^ 3; explicitating P 2 , we finally obtain 

V(P) = V K (P)+eW(P) , (3.39) 

where V K (P) = —GM/\P — 0\ is the Keplerian potential and 



h 1 

e : = 



2 4Mi? 2 



J (\Q - 0\ 2 - 3z 2 (Q)) dM(Q) . (3.41) 



In the sequel we always regard Eq. (3.39) as giving the exact potential in the region 
where the satellite is moving, and assume e > ( 2 ). The values of e when the planet 
is the Earth will be reported later on. 

By a polar motion of the satellite, we mean one with initial position and velocity 
in a plane II containing the polar axis. Such a motion stays in II for all times, since 
the gravitational force at any point of II is parallel to this plane. To analyze the 
motion we use on II a system of polar coordinates (p, •&), with ■& the angle from k. 
In these coordinates, 

C M /? 2 

W(p,0) = -==^- (1-3 cos 2 /) , (3.42) 

and the motion is described by Eqs. (3.3)-(3.5). A polar motion in the domain @ 
of Eq. (3.10) can be described via the Kepler elements (/, i?) = (P, E,Y,$), and 
we can write down Eqs. (3.22)-(3.25) or (3.33)-(3.35) for the present choice of the 
perturbation. The explicit form of Eqs. (3.33)-(3.35) is 

^ = £/*(/, [27rt]), ^ = £/«(/, [27rf]), ^ = e n/,[27rt]), (3.43) 

P(0) = P , E(0) = E , Y(0) = Y , 
depending on the function 

/ = (f )i= P ,E,Y ■■ A x T - R 3 , f(I, #) , (3.44) 



2 As an example, suppose the planet is an ellipsoid x 2 /R 2 + y 2 /R 2 + z 2 /(aR) 2 < 1 and the 
mass distribution is uniform, i.e., dM = (M/V)dV with V = (4/3)naR 3 the total volume. Then, 
passing to cylindrical coordinates r, ip, z we obtain 



AMR 2 \ J_ aR ,/n " 'I 10 



which implies e > if a < 1. As a supplementary remark, we mention that the analogous of W 
employed in [5] has a wrong sign. 



15 



/*(/, #) : = sin(^ + Y) + 2 sin(2tf) + £ sin(3tf - Y^j , 

3tt 
8P2 



37T / 

/ B (J, tf) ;= _^ 1& 2 S in(^ - 3F) + (8 + 2E 2 ) sin(tf - y) + (4 + HE 2 ) sin(tf + Y) 



+8E sin(2$ - 2Y) + 40£ sin(2$) + 2E 2 sin(3$ - 3F ) 
+ (28 + 17E 2 ) sin(3tf - F) + 24E sin(4tf - 2Y) + 5E 2 sin(5tf - 3F) ) , 

f (/, 0) := - (V cos(tf - 3F) + (8 + 6£ 2 ) cos(tf - Y) 

-(4 - 7£ 2 ) cos(tf + Y) + 8E cos(2tf - 2F) + 2AE cos(2tf) + 2£ 2 cos(3tf - 3Y) 
+ (28 + HE 2 ) cos(3tf - y) + 24£ cos(4tf - 2Y) + 5E 2 cos(5tf - 3F)) . 
This is the system to which we will apply the general scheme of Section 2. 

4 Polar motions around an oblate planet: the av- 
eraging method and its error. 

We take as a starting point the formulation of the problem in terms of the variables 
/ = (P, E, y) e A := (0, +oo) x (0, 1) x R, based on the orbit counter t as indepen- 
dent variable. So, the evolution equations have the form (2.4) with / = (f p , f E , f Y ) 
as in Eqs. (3.44). From now on, for consistency with the general notation of Sections 
1 and 2, we write 

I = (P,E,Y) (4.1) 

for the (maximal) solution of these equations for given initial data (we repeat that 
typeface symbols denote functions of t or r = et; thus I : t i— > l(t) is a function, 
to be distinguished from a point / = (P,E,Y) of the space A). Throughout the 
section, indices h, k range in {+, e, y}. 

In the sequel, we compute the averaged equations of motion and their solutions; 
then we will evaluate the error of averaging with the general method of Section 2. 

4A. The averaged system. The averages over -d G T of the functions (3.44) are 
f, where, for all / = (P, E, Y), 

T(i) = o, T(i) = o, T(i) = -%- (4.2) 

So, the averaged system for J = (J p , J E , J Y ) is 
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and has solution 
J'(r) 



3tt 



P , J E (r)=E , J y (t)=Y - forrG[0,+oo) (4.4) 

(which stays in A for all r). As we see, in the approximation given by averaging, 
the parameter and eccentricity are constant, while the argument of the pericenter 
varies linearly with the rescaled time r = si] this is a classical result. 

9f 

4B. The auxiliary functions s, v, , y, £f, Jf? . These are required by 

our general method to evaluate the error of the averaging method; their definitions 
are found in Eqs. (2.10)-(2.14). We have computed all these functions by means of 
MATHEMATICA. 

In the case of s, the components are 

s p (I, ■&) = -— 3E cos(tf + F) + 3 cos(2tf) + E cos(3tf - F) 



(4.5) 



s E {I^) 



1QP 2 



3E 2 cos(tf - 3Y) + (24 + 6E 2 ) cos(tf -Y) + (12 + 33E 2 ) cos(tf + Y) 



+ 12E cos(2tf - 2Y) + 60£ cos(2tf) + 2E 2 cos(3tf - 3Y) + (28 + 17£ 2 ) cos(3tf - Y) 
+ 18£ cos(4tf - 2Y) + 3£ 2 cos(5tf - 3F) 
1 



s y (F,tf) 



16P 2 E 



3E 2 sin(tf - 3F) + (24 + 18E 2 ) sin(tf - Y) - (12 - 21E 2 ) sin(tf + F) 



+ 12E sin(2tf - 2F) + 3QE sin(2tf) + 2£ 2 sin(3$ - 3F) + (28 + HE 2 ) sin(3tf - F) 
+ 18£sin(4tf-2F) + 3£ 2 sin(5tf - 3F) 

The expressions of v,p, q, w, u are too long to report all of them here. As examples, 
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we write down the components v p and u p , which are 
1 



v p (I, 0) 
u p (I, ■&) 



12nP 
3tt 



128P 5 L 



1QE sin Y - 18E sin(tf + Y)-9 sin(2tf) - 2E sin(3tf - Y) 
- 512P sin Y + 824P 2 sin(2F) - 103P 3 sin(tf - 3F) 



(4.6) 
(4.7) 



+ 64E 2 sin($ - 2Y) - (1296P - 36E 3 ) sin($ — F) — (768 - 1152P 2 ) sin$ 

+ (3524P + 567P 3 ) sin($ + Y) + 960P 2 sin(tf + 2F) - 24P 3 sin($ + 3Y) 

- 576P 2 sin(2tf - 2Y) + 2048P sin(2tf -Y) + (4576 + 1992P 2 ) sin(2tf) 
+ 1024P sin(2$ + Y) - 7UE 2 sin(2$ + 2Y) + 36P 3 sin(3$ - 3Y) 

+ 1216P 2 sin(3tf - 2Y) + (3180P + 521P 3 ) sin(3tf - Y) - (1792 - 640P 2 ) sin(3tf) 

- (1264P + 172P 3 ) sin(3tf + Y) - 27E 3 sin(3$ + 3Y) + 560P 2 sin(4tf - 2Y) 

- 1536P sin(4$ -Y) + (256 - 912P 2 ) sin(4$) - 336P 2 sin(4$ + 2Y) 

+ 57P 3 sin(5tf - 3Y) - 320E 2 sin(5tf - 2Y) + (288P - 144P 3 ) sin(5tf - Y) 

- (900P + 115P 3 ) sin(5$ + Y) + 88P 2 sin(6tf - 2Y) - (672 + 696P 2 ) sin(6tf) 
+ 4P 3 sin(7$ - 3Y) - (812E + 133P 3 ) sin(7$ — Y) — 328E 2 sin(8$ - 2Y) 

- 45E 3 sin(9tf - 3Y) 



For future use, it is convenient to point out that 

3tt£ 2 



P P (I) 
P E {I) 



2P 3 



■ sin(2F) , 



(4.8) 



3vr 

4P 1 

3tt 



(WE - E 3 ) sin(2F) , 



P Y (I) = ^ [34 + 25E 2 + (40 + 10E 2 ) cos(2F)] ; 



dP { ) p 3 ' dij [ ) 



otherwise ; 



^■ k (I) = for all i,j,k . 
To go on, we take & as in (2.48) and <S , M> as in (2.20); so, 

S»(l,6I,0)=Jdx^j(I + x6I,0) ] 
&j(I,6I) = [ dx^(I + xSI); 







dP 







1 QTf 

v 'dPdP 



(4.9) 
(4.10) 

(4.11) 
(4.12) 
(4.13) 
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the derivatives of s l , f l and being computed from Eqs. (4.5) (4.2) and (4. 
particular, we mention that 



In 



d 2 f 
OP 2 



(/) = 



18tt 

p4 ' 



d 2 7 

dlWl k 



(I) = otherwise. 



(4.14) 



4C. The functions R, K; time intervals. These functions are the solutions of Eqs. 
(2.21) (2.22), and can be computed in an elementary way. The expressions for R(r) 
and its inverse matrix are 



/ 1 





R(r) = 
concerning K, we find 



\ 

1 

r 1 



R-V) 



1 




\ 

1 



r(r) = 
K'(r) = 



El 
4P 

10P - 



— Tv? T 1 
p3 
-in 



67T 



(4.15) 



cos(2F ) - cos(2F --^r) 



(4.16) 



P 3 



8P 2 



67T 

cos(2Yn) - cos(2Yo - -52 r) 



K y (r) 



3tt 
167? 



34+25^+8^ cos(2F ) 



r4 



20 + P 2 



16P 2 



sin(2F ) - sin(2F - T 



In the above, r G [0, +oo); however, from now on we put a limitation r G [0, [/), for 
some finite U; consequently, the "orbit counter" t will range in [0, U/e). 

4D. The auxiliary functions /?% a 1 ^, a*., e*. fc , a 1 , 7*. These functions must 
be constructed so as to fulfill the inequalities in Section 2. First of all, we must find 
a function p = (p l ) G C([0, U), [0, +00] 3 ) such that B(j(r), p(r)) C A for all r, with 
A as in Eq. (3.27). The cited equation can be satisfied putting 



p p {r) := P , p E {r) ■= min(P , 1 - Eq) , 



p (t) := +00 



(4.17) 



for r G [0, [/) (recall that J p (t) = P and J B (r) = Eq). From these functions, we 
define T p as in Eq. (2.25); this has elements (r, r) = (r,r p ,r E ,r Y ). 
To go on, we must find a system of auxiliary functions 



a\ 0) G C°°([0,[/),[0,+oo)) , a},6\c\d},e} fc G C°°(r p , [0, +00)) 



(4.18) 



such that Eqs. (2.50) (2.51) and (2.27)-(2.30) be satisfied. As explained in subsec- 
tion 2F, each function should give an accurate upper bound for the left-hand side 
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of Eq. (2.50); for this reason we have developed an algorithm based on numerical 
maximization over t? and subsequent interpolation in r, see Appendix B. 

Let us pass to the functions a*-, b l , e* fc ; we refer again to subsection 2F, suggest- 
ing to find these functions by fairly rough majorizations of the left-hand sides of Eqs. 
(2.51) and (2.27)-(2.30). The expressions to be bounded are very lengthy in some 
cases; they were majorized with the method described in Appendix C (sometimes us- 
ing the symbolic mode of MATHEMATICA for the necessary computations). Here, 
we only report the final expressions of the majorizing functions, which are in fact 
r-independent; for this reason we write a*(r), b l (r), etc., instead of a*-(r, r), 6*(r, r), 
etc.. It should also be noted that the dependence on r is through the components 
r p , r E . In the sequel, for the sake of brevity we put 



Pn±r I 



E 



with these notations, we can take 
3 + 4P+ 



En±r' 



<{r) 



4E 4 



(4.19) 



(4.20) 



32 + 45^ + 32^ 45 + 64P+ 

— >« B ( r ) : = — ^2 — >Mn 



4 p3 



8Pi 



32 + 33P+ + 29Ej 32 + 29ffi 

a P {r)-= ,a E (r) : = o „ 2 - 2 ,a Y (r) := 



16 + 15P+ + 20^2 
32 + 30P+ + 37EI 



b p (r):= 
b E (r) : = 
b Y (r):= 
+ 

c»(r) : = 
c» : = 
+ 

c y (r): = 

+ 
+ 



4P 3 P_ 

54 + 112P+ + 33P2 

[p3 



8P1E 2 _ 



8PlE_ 



(4.21) 



6112 + 10832P+ + 6940P^ + 11372£ 3 + 1441P^ 
512PfP_ 

3520P 3 + 16384P 3 P + + 9340P 3 P^ 



256P 3 p2pf 
8940P 3 P 3 + 1861P 3 P4 

3tt 



qJ^ + T lOUli fi + t 1152P^.P 3 

504 + 1024P+ + 713P2 + 124P 3 



4608P 3 P 3 



8P 5 



(4.22) 



3tt 



148736 + 738384P+ + 1062656^ 



2048P^p2 

1220344P 3 + 675146P* + 336591^ + 26855P 
1024P 3 P 3 P 6 [ 370944P o + 2214336P 3 P + + 5434752P 3 P^ 
4927104P 3 P 3 + 2945040P 3 P^ + 1225668P 3 P^ + 147777P 3 P^ 
231936P 3 P 3 + 442368P^P 3 + 196608P^pf 
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d p Jr): = 



d' p (r):-- 



d Y p {r)-, 



e Y pp (r) 



9tt£ 2 



2P_ 



4 ' 



d p Jr) : = 



3ttE 4 



CM := ^ 



(4.23) 



p3 > 

37rP+(lO + P 2 ) 



p5 

3tt(74 + 35£ 2 



4p5 



, d Y E (r) :-- 



4 p4 
105ttP+ 

8Pi ' 



2Pi 



<(r) := 



15vr(4 + E\ 
4P 1 



18tt 



?4 ' 



e )k( T -> r ) := otherwise . 



(4.24) 



From the above functions we obtain a set of C°° functions a 1 , 7*, specializing the 
prescriptions (2.33) (2.34) to the present case. Explicitly, 



c/(r,r) := a^^r) +eb l (r) = a{ 0) (r) + <4(r)r fc + e6*(r) , 
7*(r,r,€) = 7 l (r,£) := c*(r) + dj(r)*' + £e* (r)€^ , 



(4.25) 
(4.26) 



for (r, r) G T p and £ G [0, +00) 3 . 

4E. The matrix (1 — edaz/dr) and its inverse. From Eq. (4.25), one finds 
that the derivatives of a 1 with respect to the r variables depend on r but not on r. 
From the same equation, we see that our matrix has components 



Bo 1 

5]-e^-(r)=5 



dU 



da 1 



(i,j = P , E , Y ) ; (4.27) 



all the derivatives in the right hand side are easily obtained from Eq.s (4.20) (4.21). 
Eq. (2.42), that is basic in our approach, contains the inverse of this matrix; this 
could be written explicitly, but its expression is uselessly complicated. For this 
reason, in the subsequent numerical computations we use an approximation giving 
the inverse up to a third order error Os(e, r), for (e,r) = (e,r p ,r E ,r Y ) — > ( 3 ). 
Neglecting this error, we have 

(l-e^(r)) = 1 + eM + er k N {k) + e 2 Q , (4.28) 

(with a sum for k G {p, e, y}), where we have introduced the matrices 
/ 3 + 4£ 4 4E \ 



M := 



p2 



32 + 45£ + 32£ 2 45 + 64E 16 + 15E + 20E< 



API 



8P 2 



4P 2 



32 + 33E + 29£ 2 32 + 29£ 2 32 + 30P + 37P^ 



V 



4£ P 3 



8P 2 P 2 



8P P 2 



(4.29) 



/ 



3 Rccall that, in Eq. (2.42), the matrix in which we are interested is evaluated with r % = £n*(r); 
with this position for r, (4.28) gives an inverse up to 0(e 3 ). 
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N {E) : = 



N {Y) := 



/ 



V 



4(3 + 4P 

p3 



p2 



4Po 

p2 



\ 



3(32 + 45P + 32P 2 ) 45 + 64P 16 + 15P + 20E$ 

2Pj 2P 3 2P| 

3(32 + 33P + 29P^) 32 + 33P + 58P 2 32 + 30P + 37P^ 
or d4 q zt^ d3 a tp d3 



p2 



2P 2 P 3 



45 + 64P 
2P 3 

32 + 33P °+ 58£2 



p2 




16 

p2 

16 + 29P 2 



4 

5(3 + 8P ) 



4PoPo 3 
\ 



4p2 

32 + 60P °+ 111P 2 



8P 2 P 2 



16 + 15P 



20P^ 



2P 3 

32 + 30P + 37P^ 32 + 60E + 11 IP 2 



4 

5(3 + 8P ) 
4P 2 







8P 2 P 2 



\ 



/ 



/ 



(4.30) 



Q: = 



( 746 + 1152P 



715P^ 



8P 4 



128P P 5 
64P 2 P 5 



64 + 194P + 283P 2 64 + 84P + 109P^ \ 



4P OJ P 3 



2b 



2P 3 



512P 2 P 4 



128P^P 4 



32Po-Po 
64P 2 P 4 



10208 + 27728P + 44440P 2 + 46244P 3 + 18369P 4 , 
14304 + 29344P + 71068P^ + 121568P^ + 65637P^ , 

512 + I68OP0+ ' " ' ' ' " "' 

7616 + 24832P C 



- 4405P 2 



■ 4455P^ 



3044P 4 



= 5568 + 29376P c 



25096P 2 
33400P 2 



27300P 3 
42444P 3 



7719P 4 
15153P^ 



2048 + 2880P + 7524P 2 + 5202P 3 + 4385P^ . (4.31) 

For more details on the approximate inverse (4.28), see Appendix D. 

4F. The functions R\, Pj. According to Eq. (2.35), these should fulfill the 
inequalities 



\R(t)\) < R)(r) 



R~ 



for r e [0, U) . 



(4.32) 
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On account of Eq. (4.15), we can take 



(R)(r)) := (P](r)) 



( 1 





\ 

1 



671 n 1 

T 1 



(4.33) 



/ 



4G. The main result: the estimates ^ n z (et) from Proposition 2.1. 

The "^-operation". We are finally ready to discuss the difference 

where I = (P, E, Y) is the solution of Eqs. (3.43). We follow the scheme of Proposition 
2.1; this requires to determine 4> = (^o;^o;^o) by solving a fixed point problem 



(4.35) 



In the examples that follow £ is always found numerically, as explained in Remark 
2.2(a). After £ has been found, we must solve (again numerically) a Cauchy problem 
for the unknown functions m = (m*), n = (n l ) G C 1 ([0, U), R 3 ), i.e., 



dm 1 
dr 



Ph j (en,n) , m l (0) = , 



dn l / 
rf7 = V 



(4.36) 
(4.37) 



n<(0) = £1 

with the domain conditions (2.43), having in this case the form 

< sn p < P , < en E < min(P , 1 - -Eb) , n Y > 



det 



(l- £ ^(, £ „))>0. 



(4.38) 



If the above problem has solution on [0, U), our general framework grants that the 
solution I = (P,E, Y) of (3.43) exists on [0,U/e), and gives the bounds 



\V(t)\ ^ n*(e*) for t e [0, U/e) , 



(4.39) 



i.e., 



P(t) - P | < £n p (et), I E(t) - P | ^ £n B (et), | Y(t) - J y (et)| < £tt y (et). (4.40) 
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With a slight variation of the terminology of [1] [2], we call "^-operation" the 
execution of all the numerical computations required by the present approach to 
obtain the estimators (n l ), given the initial data Pq,Eq,Eq and the values of e,U. 
These computations include: 

(a) the numerical evaluations and subsequent interpolations, necessary to determine 
the functions al s (i = p,e,y) on [0,U). We have already mentioned Appendix B, 
that describes everything in detail; 

(b) the determination of the fixed point £o; 

(c) the numerical solution of the Cauchy problem (4.36) (4.37) for (m 8 ), (n l ). 
From now on, we will indicate with %jy the CPU time required to perform 

(a)(b)(c); this time depends on the chosen hardware and software, but is anyhow 
useful for comparison with other computational approaches, such as the one de- 
scribed hereafter. 

4H. The "L-operation": a way to test the efficiency of the previous es- 
timates. This operation was already considered in [2], with a role similar to the 
one that we ascribe to it here: checking the reliability of the JV- estimates by a 
direct numerical treatment of the system (3.43), for t G [0,U/e). The direct at- 
tack of (3.43) is possible when the time scale U/e is not overwhelmingly large; for 
larger U/e, the ^-operation is too expensive, and this is just the case where the ,JV- 
procedure becomes more useful. 

The direct attack to (3.43) is performed substituting in these equations P(t) = 
P + eL p (t), E(t) = E + eL E (t), Y(t) = J Y (et) + e\/{t), which gives rise to the 
equations 

dL p 

' (t) = f p (P + eL p (t), E + eL E (t), J Y (et) + eL Y (t)\ (4.41) 



dt 
dL E 
~dT 



(t) = f E (P + eL p (t), E Q + eL E (t), J Y (et) + eL Y (t)), 



L> -(/! = f Y (P + eL p (t), E + eL E (t), J Y (et) + eL Y (t)) + ? " n 



dt y j J v u w ' u v y ' v 7 y " J p (ety ' 

(L-,L B ,L y )(0) = 0. 

By definition, the L-operation is the numerical solution of this Cauchy problem in 
the unknowns V, for t G [0, U/e); will indicate the necessary CPU time. 

5 Examples. 

Introducing the examples. In the case of the Earth, 

m 

GM = 3.98600442 • 10 14 — -, R = 6.378135 • 10 6 m , (5.1) 

sec^ 
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e = 5.457 • 10" 4 



So, for a satellite on a polar orbit around the Earth, the unperturbed apogee and 
perigee p + , p_ and the orbital period T orb are determined in this way by the initial 
data (P ,Eo,Y ): 

p p3/2 

p T = Q — x 6.378135 x 10 3 Km, T orb = - ^— x 1.408150 hours. 

1 ± E (1 — EqY' 2 

In the sequel we present two examples, with initial data (Pq, Eq, Yq) very similar 
to the actual data of two real satellites (Polar and Cos-B). The MATHEMATICA 
package, already mentioned in relation to the symbolic treatment of the problem, 
has been employed (on a PC) for the required numerical computations. 
We repeat a comment of the Introduction: the examples are not fully realistic 
since they do not account for perturbations different from the J2 gravitational term; 
nevertheless, we think they have some interest because they show the effectiveness 
of the method, suggesting that our approach could be applied as well including other 
perturbations. 

Each example is worked out along these lines: 

(i) Firstly, we take U/e = 3000. We perform both the JV- and the -^-operations, 
using the second to test the first one. We give figures reporting the graphs of e\ V(t)\ 
(rapidly oscillating) and eri l (et), for i E {p,e,y} and t E [0, U/e); the CPU times 

%c are also indicated. 

(ii) Next, we choose U/e = 60000. The ^/T-operation is still performed within short 
CPU times; we give figures reporting %jy and the graphs of the estimators (en*). 
On the contrary, the ^-operation exceeds the capabilities of the machine employed. 

Example 1. We take 

P := 3.000 , E := 0.6640 , Y := 0.0000 , (5.2) 
corresponding to 

p+ = 56950 Km , p_ = 11500 Km , T orb = 17.50 hours . (5.3) 

These data are very similar to the initial conditions of the Polar satellite, taking 
t = on April 1997 [8]. Eq. (4.4) for the solution of the averaged system on any 
interval [0, U) gives J p (t) = const. = Po, J b (t) = const. = E Q , and 

j Y (r) = -1.047r . (5.4) 

Figures la, lb and lc report the graphs of n l (et) and \V(t)\ for t E [0,3000], as 
produced by the £- and ^K-operations. We note that 3000 T orb ~ 6 years. The CPU 
times (in seconds) for jV and £ reported in Fig. la refer to the overall calculation, 
also including the e and y components plotted in Figures lb, lc. 
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Figures Id, le and If report the graphs of vC(et) and | V(t)\ for t G [0,60000], as 
produced by the ^F-operations. We note that 60000 T orb ~ 120 years. The CPU 
time indicated in Fig. Id also includes calculations of the e and y components, plotted 
in Figures le, If. In comparison with the previous case t G [0,3000], the number of 
orbits is increased by a factor 20; in spite of this, T jy changes by a factor < 2. 

Example 2. We take 

P := 1-973 , E = 0.8817 , Y := 0.9600 , (5.5) 

corresponding to 

p+ = 106400 Km , p_ = 6688 Km , T orb = 37.16 hours . (5.6) 

These data are very similar to the initial conditions at launch (August 1975) of the 
Cos-B satellite [9]. Eq. (4.4) for the solution of the averaged system on any interval 
[0, U) gives J p (r) = const. = P , J E (r) = const. = E and 

j Y (r) = 0.9600 - 2.421 r . (5.7) 

All figures for this example give the same information as the corresponding ones 
of Example 1. In particular, Figures 2a, 2b and 2c refer to computations for t G 
[0,3000], while Figures 2d, 2e, 2f refer to the case t G [0,60000]. We note that 
3000 T orb ~ 12 years and 60000 T orb ~ 250 years. The CPU times (in seconds) for 
,jV and £ are similar to the ones in Example 1, so the comments given therein also 
apply to this case. 
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Figure la. Example 1, for t E [0,3000]. 
Total CPU times: 7> = 5.18 sec, T L = 
29.08 sec. Graphs of en p (et) and e\L p (t)\. 
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Figure Id. Example 1 , for t G [0, 60000]. 
Total CPU time: 7> = 10.11 sec. Graph 
of en p (et). 
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Figure lb. The same situation as in Fig. Figure le. The same situation as in Fig. 



la. Graphs of en E (et) and e\L E (t)\. 



Id. Graph of en E (et). 
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Figure lc. The same situation as in Fig. Figure If. The same situation as in Fig. 



la. Graphs of en Y (et) and e|L y (t)|. 



Id. Graph of en Y (et). (Note that J Y (et) 
-34.29 for t = 60000). 
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Figure 2a. Example 2, for t £ [0,3000]. 
Total CPU times: 7> = 5.77 sec, T L = 
33.02 sec. Graphs of en p (et) and e\L p (t)\. 
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Figure 2d. Example 2, for t £ [0, 60000]. 
Total CPU time: 7> = 11.28 sec. Graph 
of en p (et). 
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Figure 2b. The same situation as in Fig. Figure 2e. The same situation as in Fig. 



2a. Graphs of en E (et) and e\L E (t)\. 



2d. Graph of en E (et). 
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Figure 2c. The same situation as in Fig. Figure 2f. The same situation as in Fig. 



2a. Graphs of en Y (et) and e\L Y (et)\. 



2d. Graph of en Y {st). (Note that J Y (et) 
-78.31 for t = 60000) 
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A Appendix. On Kepler elements. 



We refer to the framework of subsection 3B, having fixed = (p, p, i?) G f^. 

The unperturbed Kepler orbit with these initial data at time t = is an ellipse 
with (dimensionless) parameter P = P(£, i?), eccentricity P = P(£, $) and argument 
of the pericenter Y = Y(£, #). In any textbook on classical mechanics, one finds the 
equations 

— ( \ -\\ 

P ~ l + Pcos(tf- Y) ' 1 j 

A pppsin^-Y)^ 



;i + Pcos(^-r)) 2 



Eq. (A.l) is the representation of a conic in polar coordinates; Eq. (A. 2) follows 
taking the derivative of this polar representation with respect to time along the 
unperturbed Kepler motion, i.e., regarding P, E, Y as time independent. Eqs. (A. 3) 
(A.4) give the parameter and the eccentricity as functions of the angular momentum 
and the total energy (per unit mass) p 2 $ and (l/2)p 2 + (1/2) p 2 $ 2 — GM/p . 

We note that Eq. (A.l) and Eqs. (A. 3) (A.4) are equivalent to Eq. (A.l), to the 
square of Eq.(A.2) and to Eq. (A. 3). Eq. (A. 2) is important because (due to d > 
in Qi) it implies 

signp = signsin(-# — Y) ; (A. 5) 

this information is essential to determine Y as a function of (£, i?). Eqs. (A.1)-(A.5) 
yield the expressions for P, E, Y in Eqs. (3.11)-(3.13); from here, one also proves the 
bijectivity of the map (£, i?) i— > (P, E, Y, i?) = (/, i?) and the expression (3.16)-(3.19) 
for its inverse. 

To conclude, let us comment on the equations (3.20) for the apocenter p + , the 
pericenter p_ and the orbital period T or b of the unperturbed Kepler motion along an 
ellipse of parameter P and eccentricity E. The expressions for p± are consequences 
of the polar equation for the ellipse; the expression for T orb follows from Kepler's 
third law 

rp 2tt . . 3/2 2tt (p+ p_\3/2 

i or b = . (semi-maior axis) 7 = . (A.o) 

VGM V ; VGM V 2 2 J K J 

and from the previously mentioned formulas for p±. 
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B Appendix. Numerical computation of the func- 



tions a 



(0)' 



Let us consider a general periodic system (2.4), and suppose we have the related 
functions s, J, R, K, for r within a finite interval [0, U). Hereafter we outline a simple, 
but effective scheme to construct by mixed, numerical and analytical, techniques a 
function r G [0,U) i— > a l , Jr) fulfilling up to small errors the inequality (2.50), for 
any given % G {l,...,d}; this method is easily implemented on MATHEMATICA. 
The scheme consists of four items, 
(i) Take on the torus T a grid of equally spaced angles 



V= [^r] > (g = i,2,...,Q) . 



fB.l^ 



r„:=-n 



(ii) Take in [0, U) a grid of equally spaced instants 

n = 0,1,..., AT- 1 . 

(iii) For each n, compute numerically 

;j(r n ), # q ) - R(r n )s(/ , # ) - K(r n ))' 



(B.2) 



(B.3) 



then, for all d G T, 



s(j(r„),tf) -R(t)s(J ,#o) -K(r n ) 



< a 



(0)n 



(B.4) 



up to an error neglected in the sequel, which is of order 0(1 /Q 2 ) for Q — > +oo. 
(iv) Interpolate the sequence a % n (n = 0,...,N), finding a smooth function r G 
[0, U) I— > a( )( r ) sucn 



°(o) ( r «) = a 



(0)n 



for n = 0,1,..., N-1. 



(B.5) 



More precisely, we define d 1 ,^ to be the Lagrange polynomial (restricted to [0, U)) 
such that O( ) (r n ) = O( ) n for all n; thus 



AT-l 



AT-l 



a (0) (r 



) = E a (0)n ( II f^T") for TG [0,f/), 



(B.6) 



n=0 



m = 
m ^ n 
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with r l n and a l ^ n as in Eqs. (B.2) (B.3). By construction, the inequality (2.50) 

(s(J(r),7?)-R(r)s(/ ,0)-K(r)) 

is fulfilled, up to a small error, at all points r = r n ; we assume the same to happen 
everywhere in [0, U). Of course, being C°°, the functions (B.6) fulfill the regularity 
conditions for application of the basic Proposition 2.1. 

The above scheme (i)...(iv) has been employed to construct the functions a 1 ^ 
(i = p, e, y) in the numerical examples of Section 5 on the motions of satellites. In 
these examples, Q = 30 and N = 100 (with such an N, the Lagrange interpolating 
polynomials are computed very quickly by MATHEMATICA). Admittedly, the ap- 
proach followed does not account for the small errors mentioned in items (iii) (iv): 
throughout the paper these errors, and the ones produced by the approximate in- 
version of the matrix 1 — eda/dr, are the only ones for which an analytical estimate 
is not provided. 



<a{ 0) (r) for all?? GT 



C Appendix. The functions a*, 6%..., e l - k for the 
satellite problem. 

As usually, the indices k range in {p, e, y}. 

Finding a*.. As an example, let us illustrate the determination of a p . To this 
purpose, we use Eq. (4.11) with i = p,j — y, I — J(r) = (P , E , J Y (r)) and 
51 — 5 J, giving 

f 1 ds p 

S?(J(t),6J,0)= dx—(P + x5J p ,E + x5J E ,J Y {T)+x5J Y ,$) (C.l) 
Jo 

1 p _|_ T A T E / \ 

dx (3sm($ + J Y (T)+x5J Y )-smm-J Y (T)-x5J Y )) ; 

P + x5J p V / 

in the last passage, the derivative ds p /BY has been computed from Eq. (4.5). Now 
we take the absolute value in the last equation, and use the elementary inequalities 
\^dx...\ <: fodx\...\, \E + x5J E \ ^E + \5J E \, l/\P + x5J p \ < 1/(P - |^ P |), 
| sin(..)| ^ 1. In this way we get 

\y p (j(r),SJ^)\ < 4 ^ + 1^1 = a p (r)\ ri=lsjil , (C.2) 

where a p (r) := 4(£o + 5r E )/(Po — r p ); this definition of a p agrees with Eq. (4.20) 
(where we use the abbreviations E + := E + r E , P_ := P — r p ). 
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Finding b x and c\ We illustrate the determination of b Y . To this purpose, we 
must first compute the p component of the function on the left-hand side of Eq. 
(2.27) (for J = (J*) and 5 J = (5 J 1 ) (i = p,e,y). Recalling again that J p (r) = P 
and J b (t) = E Q for all r, one finds 



{w( J(r) + 5 J, 0) - ^( J(r)) „( J(r) + 5 J, 0))' 



(C.3) 



= T7 777 77 P^ sin(A;tf + £(j y (r) +cLT)) , 

5i2P 3 (p 0+ ^r(p 0+ ^ ^ ^_ 6 _ 6 ke 

where the coefficients Pu are polynomials in P , P + <5J P , P + SJ E ] for example, 

P 01 = 256(P + 5J E )(-49P 3 + 13P 3 (P + 5J B ) 2 -16(P + 5J P ) 3 (Po + 5J E ) 2 ) . (C.4) 
This implies 

97, 



(w(J(t) + 5 J, 0) - ^( J(r)) V ( J(r) + 5 J, 0)) 



E 



(C.5) 



512(P -|5J p |) 3 (Po-|^ B |) 2 



fc=0,...,10;^=-6,...,6 

the absolute value of each coefficient Pu is easily bounded, using fairly rough esti- 
mates such as |P i| < 256(Po+|5^|)(49P 3 +13P 3 (P +|5^|) 2 +16(Po+|5J P |) 3 (Po+ 
\5J E \) 2 ). The conclusion is 



w(J(r) + 5 J, 0) - ^ ( J(r)) W ( J(r) + 5J, 0) 



^ b Y (r) 



|<5J l 



(C.6) 



where 6 Y is the function appearing in Eq. (4.21) (which is independent of r Y ). The 
determination of b p ,b E and c p ,c E ,c Y is performed similarly, yielding Eqs. (4.21) 
(4.22). 

Finding d x y As an example, let us consider the determination of d p . We use Eq. 
(4.12) with i = p,j = Y ,I = J(r) = (P , E , J y (r)) and 51 = 5 J, giving 

'(J(t), <JJ) = [dx^(E + x5J E , P + x5J p , J y (t) + x5J y ) (C.7) 
Jo ()Y 

"\ {E + x5J E ^ 2 



— 37T / dx 



in the last passage, the derivative &p p /dY has been computed from Eq. (4.8). The 
last equation implies 



(E +\5J E \) 2 

(P -\5J p \f "^''irHw*! 
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J(t),5J)\ ^ 3tt 



dUr 



(C.8) 



with d p Y as in Eq. (4.23). 



Finding e*. fc . These functions should bind the absolute values of the components 
J^4(j(r), 8 J) of JT; these are provided by Eq. (4.13), with / = J(r) and 51 = 5 J. 
From the cited equation, we see that we can take e* fc (r, r) := for (i, j, k) ^ {y, p, p). 
Again from (4.13), we infer 



^ 36tt 



je p Y p (j(T),5J)\ = 36tt 



f 

Jo 



dx 



:i-x) 



(P + x5J P Y 



x 



18tt 



(Po-|^1) 4 {Po~\5J p \Y 
with e pp as in Eq. (4.24) . 



e Y PP {r) 



r i =\5J i 



(C.9) 



D Appendix. The approximate inverse (4.28). 

We must invert the 3x3 matrix 1 — e(da/dr)(r), whose elements are given by Eq. 
(4.27). For (e,r) — > 0, this equation implies 

l-e-£(r) = l-eM -er k N {k) -e 2 Z + 3 (e,r) for(e,r)->0, (D.l) 

where we provisionally write M, and Z for the matrices of elements 

(9(7* da 1 r)U 

M):=a}(0), Ar (fc) ;.:=^(0) + ^(0), Z) := — (0) , (D.2) 

(i,j,k G {p,e,y}). Eq. (D.l), with the standard X — > expansion (1 — X) -1 = 
1 + X + X 2 + 3 (X), implies (1 - e{da / dr)^))' 1 = 1 + (eM + er k N (k) + e 2 Z) 
+ (eM + er k N {k) + e 2 Z) 2 + 3 (e, r) = 1 + eM + er k N (k) + e 2 Z + e 2 M 2 + 3 (e, r). 
Summing up, 

(l-e^(r)) = l + eM + er k N {k) + e 2 Q + 3 (e,r) , Q := Z + M 2 . (D.3) 

This justifies the approximate inverse (4.28) if we check that the matrices M, 
and Q, defined here via Eq.s (D.2) (D.3), coincide with the homonymous matrices 
in Eq.s (4.29) (4.30) (4.31). This is obtained substituting in (D.2) (D.3) the explicit 
expressions (4.20) (4.21) of the functions a* and b l . 
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